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Ultralong-range order in the Fermi-Hubbard model with long-range interactions 
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We use the dual boson approach to reveal the phase diagram of the Fermi-Hubbard model with 
long-range dipole-dipole interactions. By using a large-scale finite-temperature calculation on a 
64 X 64 square lattice we demonstrate the existence of a novel phase, possessing an ‘ultralong- 
range’ order. The fingerprint of this phase - the density correlation function - features a non-trivial 
behavior on a scale of tens of the lattice sites. We study the properties and the stability of the 
ultralong-range ordered phase, and show that it is accessible in modern experiments with ultracold 
polar molecules and magnetic atoms. 
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Recent experimental progress opened up a possibil¬ 
ity to use ultracold quantum gases in optical lattices 
to realise exotic many-particle Hamiltonians inaccessi¬ 
ble in ‘conventional’ condensed matter physics^^^. One 
of the rapidly advancing research directions deals with 
particles possessing a dipole moment, such as magnetic 
atoms^^® or ground-state heteronuclear molecules 
The anisotropic and long-range character of the dipole- 
dipole interactions between such species is predicted to 
give rise to novel phases of matter^’To date, sev¬ 
eral many-body models have already been implemented 
in the laboratory^^’^^. 

In particular, optical lattice experiments allow to sim¬ 
ulate the Dipolar Fermi-Hubbard (DFH) model, i.e. an 
extended Hubbard model with long-range dipole-dipole 
interactions^^. While the phase diagram of its bosonic 
counterpart has been evaluated using large-scale quan¬ 
tum Monte Carlo simulations^®’^^, understanding the 
DFH model represents a formidable challenge due to the 
sign problem in quantum Monte Carlo^®’^®. 

Recently, the DFH model has been approached us¬ 
ing a number of techniques, starting from traditional 
mean-held^°“^^ and Hartree-Fock-Bogoliubov approxi¬ 
mations^^. The Functional Renormalization Group tech¬ 
nique^^ has uncovered novel bond-ordered phases in sys¬ 
tems of dipolar fermions at half-filling^^’^®. 

However, it is unclear how applicable it is when the 
local and dipole-dipole interaction strengths are compa¬ 
rable with the kinetic energy. Furthermore, most optical 
lattice experiments take place in a harmonic trap, with 
different kinds of hllings present simultaneously. There¬ 
fore, it is important to understand the behavior of the 
system away from the special case of half-filling. 

Dynamical Mean Field Theory (DMFT)^^’^® has been 
extensively used for electronic structure theory^®. It can 
be applied both at small and at large local interaction 
strength. This method has been adapted to the study 
of optical lattices in the form of real-space DMFT^°’^^ 
and used to study phase transitions in half-filled dipo¬ 
lar fermion systems^^. However, the nonlocal interac¬ 
tion had to be restricted to its Hartree contribution. In 
real-space DMFT a single-site problem has to be solved 
for every lattice site, hindering the study of large lattice 


sizes. 

Here, we reveal the phase diagram of the DFH model 
using the dual boson approach^^. The method is based 
on a single-site impurity problem and therefore allows to 
treat lattices of larger sizes, which is crucial for systems 
featuring long-range interactions. On the other hand, 
the technique is applicable even when the interaction 
strength is comparable to the kinetic energy. Further¬ 
more, the technique allows to perform finite-temperature 
calculations away from half-filling, which is crucial in or¬ 
der to reproduce the conditions of realistic experiments. 
As the main result, our large-scale calculation allows to 
demonstrate the occurrence of a novel phase, featuring 
‘ultralong-range’ density correlations at distances of tens 
of lattice sites. While such a phase has never been pre¬ 
dicted before, it appears to be within the reach in modern 
experiments with dipolar quantum gases. 

We start with a two-component gas of dipolar fermions 
trapped in a two-dimensional square optical lattice, as 
schematically illustrated in Fig. 1, and described by the 
DFH Hamiltonian: 

H = Y X ^ X (1) 

(jk)a j jk 

Here creation (annihilation) operator for 

a fermion with spin state a on site j, and {jk) is a pair 



Figure 1. (Color online) Dipolar fermions in a square opti¬ 
cal lattice. The orientation of the dipoles is given in by the 
spherical angles {0,(j)). 
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of nearest neighbours. In experiment, the two spin com¬ 
ponents cr can be represented by different fine, hyper- 
fine, or rotational states of the dipolar species used. The 
first two terms of eq. (1) give the amplitudes of the 
nearest-neighbor hopping, t, and the on-site interaction, 
U. The third term corresponds to the dipole-dipole in¬ 
teraction, whose spatial dependence is given by = 


Cd 1 — 3(fjfc ■ /i^jk/a)^- Here rjk is the vector con¬ 

necting the fermions on sites j and fc, with fjk = rjk/rjk, 
and a is the lattice constant. The fermions’ dipole mo¬ 
ments point in the same direction given by the unit vec¬ 
tor d. The dipole-dipole interaction strength parameter 
depends on the particular species involved and is given 
by Cd = d^/(47reo) for the electric dipoles of magnitude 
d in the laboratory frame, and by Cd = /(in) for 

the magnetic dipoles of magnitude /x. Here Cq and 
give, respectively, the permittivity and permeability of 
vacuum. 


Table I illustrates the strength of the dipole-dipole 
interaction parameters that can be achieved with sev¬ 
eral ultracold species currently available in the labo¬ 
ratory. For the molecules, the dipoles can be conve¬ 
niently polarized using a microwave held coupling the 
two lowest rotational states, J = 0 and J = 124,44^ 
this case, due to the contributing Clebsch-Gordan co¬ 
efficients, the resulting magnitude of the dipole-dipole 
interactions of Table I needs to be divided by a factor 
of six. For a static held with a magnitude E, the ra¬ 
tio, d/dmoi, of the induced dipole moment to the molec¬ 
ular one can be estimated in the strong-held limit as 
d/dmoi = [l ~ {2dmoiE/B)~^/‘^Y^, where B is the molec¬ 
ular rotational constant. Experimentally feasible helds 
thus allow to achieve d/d^oi ~ 0.7 — 0.8^®, which results 
in the reduction of the dipole-dipole interaction matrix 
elements by a factor of two compared to the values listed 
in Table I. We see that for typical values of lattice hop¬ 
pings, t ~ 10 — 10® Hz, values of Cd > t are achievable for 
dipolar molecules; furthermore, the regime of substantial 
magnitudes of Cd can be accessible with magnetic atoms. 
We note, moreover, that recently created ultracold Er 2 
molecules^®, if prepared in their fermionic incarnation, 
can in principle allow for Cd = 135.2 Hz on a lattice 
with a = 266 nm. In experiment, the orientation of the 
dipoles with respect to the lattice, as given by the angles 
{6, (j)) of Fig. 1, can be controlled by tilting the polariza¬ 
tion of the external microwave, electrostatic, or magnetic 
held®’i^-i®. 

The main idea of the dual boson method^® is to sepa¬ 
rate the local and nonlocal physics from each other. As 
in DMFT, the local physics is encapsulated in an aux¬ 
ilary single-site impurity problem. Since the impurity 
problem possesses only a few degrees of freedom, it can 
be solved numerically exactly even when the correlation 
effects are significant. The remaining nonlocal physics 
is decoupled using an exact transformation to the new - 
dual - degrees of freedom. The main correlation effects 
are taken into account at the level of the impurity prob- 


Table I. The dipole-dipole interaction strength, Cd (in Hz), 
for selected fermionic species currently available in the labo¬ 
ratory. For molecules, the value of the molecular-frame dipole 
moment was used to evaluate Cd- 
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lem, thus the dual degrees of freedom are only weakly cor¬ 
related, and therefore can be treated using perturbation 
theory. The lattice size enters only the relatively simple 
dual part of the calculation, which allows to increase the 
system size at a reasonable computational cost. In turn, 
this allows to signihcantly increase the momentum reso¬ 
lution, which is required for studying long-range ordered 
phases. As an example, the method has recently been ap¬ 
plied to plasmons in two-dimensional strongly-correlated 
electron systems^^, where the long-range Coulomb inter¬ 
action plays a crucial role"*^^. More details about the com¬ 
putation scheme can be found in Refs.^®’^®. 

In the dual boson approach, a charge order instability 
is revealed as a divergence in the static density-density 
susceptibility, Xq = (pp)q^=o> p = n — (n). The 
latter can be obtained from the dual perturbation theory 
using an exact transformation^®, and satisfies the charge- 
conservation laws®®. The momentum q, in turn, charac¬ 
terizes the type of the emerging charge order. This way, 
the signatures of a phase transition are already visible as 
it is approached from an unordered phase. In experiment, 
this corresponds to real-space density correlations which 
emerge in the vicinity of the phase transition. In an ul¬ 
tracold setting, the density-density correlations can be 
detected using, e.g., Bragg scattering®®, time-of- flight®^, 
or noise correlation®®’®^ spectroscopy. 

We now evaluate the phase diagram of the DFH Hamil¬ 
tonian (1) with a repulsive on-site interaction U = 4t, the 
dipole-dipole interaction strength Cd = 2t, and filling of 
(n) « 0.9 fermions per site, with equal populations of 
spin-up and spin-down states. In order to reproduce the 
conditions of ultracold experiments we set the tempera¬ 
ture to T = t/{AkB), where fcs is Boltzmann’s constant. 
As an example, for the hopping rate of t = 2TTh x I kHz 
this corresponds to T = 12.5 nK. In order to reveal the ef¬ 
fect of anisotropy on the many-body state of the system, 
we evaluate the phase diagram depending on the orienta¬ 
tion of the dipoles with respect to the lattice plane. The 
resulting phase diagram is shown in Fig. 2. 

Depending on the orientation of the dipoles, the sys¬ 
tem goes through different phases. In the middle of the 
phase diagram, there is a normal phase featuring no par¬ 
ticular order (a “metal” or “Fermi liquid” state). At 
smaller values of 0, which corresponds to the dipoles ori¬ 
ented nearly perpendicular to the lattice plane, a tran- 
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Figure 2. (Color online) Phase diagram as a function of the 
dipole orientation at 17 = 4f, Cd = 2t, and filling (n) « 0.9. 
The dotted lines show the phase boundaries at the reduced 
dipolar coupling, Cd = 1.8t. The black diamonds show the 
angles selected for Fig. 4. 
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Figure 3. (Color online) Inverse charge susceptibility in the 
normal phase, for the same parameters as in Fig. 2. The 
blue triangles show the divergence of the susceptibility at the 
checkerboard point, qcs = (7r,7r), as Q is lowered, both for 
(j> = 0 (filled triangles) and (j) = 0.257r (empty triangles). 
The green squares show the divergence of the qstripe = (0, tt) 
susceptibility as 6 increases ai (f) = 0 and the red pentagons 
show the q, « (0.2l7r, —0.2l7r) susceptibility at <)> = 0.257r. 
The dashed lines show a linear extrapolation of the inverse 
susceptibility. 
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Figure 4. (Color online) (a) Momentum-space susceptibility 
at selected points of the phase diagram. Fig. 2. (b) The cor¬ 
responding density correlation function in real-space: given a 
particle in the center of the figure, red indicates a higher prob¬ 
ability to find a particle at x, y and blue a lower probability. 
Every pixel corresponds to a lattice site. 
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Figure 5. (Color online) Real space density correlation func¬ 
tion as in Fig. 4 (b), at fixed 9 = O.SStt, for three values of cj>. 
The areas of high and low density follow the angle 4> (black 
line). 


sition to checkerboard order [with qcB = (tTiTt)] oc¬ 
curs. Such a checkerboard phase has been previously 
observed in half-filled systems with nearest-neighbor in- 
teraction'^®’®®“®° and away from half-filling at 0 = 0^^. In 
the normal phase, the signatures of the transition are al¬ 
ready visible as the phase boundary of the checkerboard 
phase is approached: as shown in Fig. 3, the suscep¬ 
tibility at qcB = (tt, tt) diverges as 9 is lowered. The 
susceptibility is shown both at ^ = 0 (filled symbols) 
and (j) = 0.257r (empty symbols), the checkerboard sus¬ 
ceptibility depends on (j) only weakly. The left panels of 
Fig. 4 show the charge susceptibility close to the phase 
boundary, at 0 = 0.127r, (j) = Q. The momentum-space 
susceptibility has a clear maximum at qcB = (tt, tt). The 
sign of the real-space charge correlation function, shown 
in panel (b), features a checkerboard pattern. 

At larger 9, when dipoles are oriented nearly in the lat¬ 
tice plane, other charge-ordered phases occur in Fig. 2. 
At small (j) (dipoles pointing along x), there is a horizon¬ 


tally striped phase, whose charge order is given by the 
momentum qstripe = (0,7r). Again, the tendency towards 
a diverging susceptibility is already visible in the normal 
phase, see Fig. 3. The real and momentum space suscep¬ 
tibility close to the phase boundary, at 9 — 0.207r, (j) = 0, 
shows a stripe pattern at small x. As the phase bound¬ 
ary is approached, the striped order takes over also at 
longer wavelengths (see Supplemental Material). A sim¬ 
ilar striped phase has been predicted before^®. 

Finally, when both 9 and (j) are large, corresponding 
to the dipoles oriented along the aiy-diagonal, we find 
a novel phase possessing an ultralong-range order. The 
right-hand side of Fig. 4 shows the susceptibility close 
to the transition towards this ordered phase. Contin¬ 
uous areas of high and low density extend over a large 
number of lattice sites. The maximum of the momentum- 
space susceptibility also occurs at smaller q (longer wave¬ 
length). Unlike in the other two ordered phases, here the 
shape of the real-space correlation function strongly de- 
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pends on the angle with the areas of high and low 
density rotating along with the dipole orientation in the 
plane. The ordering vector q* depends on (p and evolves 
from q* « qstripe = (0, —tt) close to the border with 
the striped phase to q* « (0.2l7r, — 0.2l7r) at (/> = 0.257r. 
Fig. 5 shows that the high density parts of the real-space 
susceptibility follow the dipole angle 4>. The angular de¬ 
pendence of the density correlation function reflects the 
anisotropy of the dipole-dipole interaction, as given by 
the second spherical harmonic. 

The spatial symmetry of the system helps to get in¬ 
sight into the different orderings. When the dipoles are 
oriented in the z-direction [0 = 0], the system possesses 
mirror symmetry with respect to the x- and y-axis and 
the xj/-diagonal, which coincides with the symmetry of 
the checkerboard pattern. For dipoles oriented along the 
x-axis [9 = 7r/2, p = 0], the mirror symmetry along 
the diagonal is broken by the dipole orientation. Indeed, 
the horizontally-striped phase has the mirror symme¬ 
tries only with respect to the x- and y-axes. Finally, for 
dipoles oriented diagonally in-plane [9 = 7r/2, p = 7r/4], 
the system has mirror symmetry around the two in¬ 
plane diagonals - exactly as the resulting ultralong-range 
charge order. 

Let us discuss the stability of the observed phase di¬ 
agram with respect to changing the parameters of the 
Hamiltonian (1). The ordered phases originate from the 
interplay between the contact and the long-range interac¬ 
tion terms, therefore the phase boundaries shift depend¬ 
ing on the value of c^. This point is illustrated in Fig. 2 
by the dotted lines, corresponding to the phase bound¬ 
aries at a reduced dipolar interaction, Cd = 1.8t. The 
structure of the phase diagram, however, stays qualita¬ 
tively similar. In addition, we have studied the phase 
diagram as a function of the filling (n). At (n) « 0.8, the 
phase diagram is qualitatively similar to Fig. 2, while at 
a further reduced value of (n) « 0.49 the checkerboard 
instability disappears at small 9. In Fig. 2, the normal, 
striped and ultralong-range phases meet near (p = 0.127r. 
The location of the triple point depends on the density 
(n) and the dipole strength c^, however we found it al¬ 
ways to be close to ^ = 0.127r. This indicates that in 
experiment small density fluctuations due to a harmonic 
trap are unlikely to qualitatively change the observed or¬ 
der. Furthermore, ultracold gases in uniform potentials 
with constant density have been created in laboratory®^. 

In order to study the importance of the long-range 
character of the dipole-dipole interactions, we have per¬ 
formed simulations with the interaction 14“^ restricted 
to the nearest neighbours. In such a case, the suscep¬ 
tibility near the boundaries with the checkerboard and 
striped phases looks qualitatively similar. On the other 
hand, the dumbbell-shaped density correlation function 
that occurred in Fig. 4 at 0 = O.SStt, ip = 0.257r, dis¬ 
appears completely. Therefore, we conclude that while 
the checkerboard and striped phases are mainly driven 
by the nearest-neighbor interaction, the long-range cou¬ 
plings are essential for the ultralong-range ordered phase 


to form. 

Large-scale simulations appear to be vital in order to 
capture and describe the ultralong-range order. As il¬ 
lustrated in the Supplemental Material^®, in simulations 
involving smaller lattice sizes finite-size effects can signif¬ 
icantly affect the results. 

Thus, we have demonstrated the occurrence of a novel, 
ultralong-range-ordered phase in the Fermi-Hubbard 
model with dipole-dipole interactions. This reveals 
the importance of taking into account large system 
sizes while dealing with the species featuring anisotropic 
long-range interactions. The novel phase should be 
within reach for current experiments with ultracold po¬ 
lar molecules and magnetic atoms. The formation of the 
ultralong-range order might be a general phenomenon 
for systems with interactions beyond the nearest neigh¬ 
bor, e.g. those involving quadrupoles®^^®^ and oscillating 
light-induced dipoles®®’®®. 
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Appendix A: Finite-size effects 

The simulations were performed on a 64 x 64 square lat¬ 
tice with periodic boundary conditions. Selected points 
in the phase diagram were also calculated on a 128 x 128 
lattice in order to rule out the finite-size effects. For 
comparison, we have also done calculations on a 16 x 16 
lattice. There, finite size effects are clearly visible and the 
momentum resolution is insufficient to accurately deter¬ 
mine the ultralong-range order. The comparison of the 
results for different lattice sizes is shown in Fig. 6. 


Appendix B: Density 

The dual boson approach works in the grand canoni¬ 
cal ensemble, therefore all results are obtained at a hxed 
chemical potential. We fixed the chemical potential af¬ 
ter subtracting the Hartree contribution of the nonlo¬ 
cal interaction. We note that while scanning over the 
dipole orientation, as in Fig. 2 of the main text, the den¬ 
sity cannot be held exactly constant since the dipolar in¬ 
teraction also gives beyond-Hartree contributions to the 
density. However, the difference in density between dif¬ 
ferent points in the phase diagram is on the order of 
A(n) « 0.01 and is therefore negligible. 
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Figure 6. (Color online) The real-space correlation function at 
6 = O.SStt, ij> = 0.25n, for three different lattice sizes. Values 
are rounded to the range of [—0.01,0.01]. The 64x64 and 
128x128 results are very similar; the 16x16 result, on the 
other hand, reveals finite size effects. 


Appendix C: Computational scheme 

In the calculation of the susceptibility, we first calcu¬ 
lated the susceptibility of the dual particles by summing 
all particle-hole ladder diagrams. In this way, we treated 
repeated scattering processes to all orders and the re¬ 
sulting susceptibility satishes the charge-conservation re- 
quirements^^’^®. 

Charge-order transitions can manifest themselves in 


two ways in our calculations. As explained in the main 
text, a divergence in the dual boson charge suscepti¬ 
bility is a sign of an ordered phase. However, some¬ 
times the EDMFT self-consistency (which is done prior 
to the dual boson part) already shows a diverging sus¬ 
ceptibility. In such a case, converging to EDMFT self- 
consistently is difficult. Previous work on the checker¬ 
board charge-ordering transitions in the extended Hub¬ 
bard model showed that the dual boson transition occurs 
already at a smaller interaction strength compared to the 
EDMFT transition^®. However, the region between the 
dual boson ordering and EDMFT ordering can be small. 
For the results of Fig. 2 of the main text, EDMFT con¬ 
verged in the checkerboard region, and the phase transi¬ 
tion became visible when the nonlocal dual boson correc¬ 
tions were added to the susceptibility. In the striped and 
ultralong-range ordered phase at large 9, on the other 
hand, the EDMFT part of the calculation already suffers 
from divergencies. Here, following the susceptibility from 
the normal phase is the most stable way to determine 
the type of order. The location of the triple point of the 
normal phase, where it meets with the striped and ultra¬ 
long-range order, can be determined by approaching it 
from the normal phase. However, using our method, it 
is challenging to predict how the the phase boundary 
between the striped and ultralong-range ordered phases 
behaves for angles 6 above the triple point. 

To study the properties of the normal phase (with¬ 
out order), we have looked at the single-particle Green’s 
function in imaginary time, G(t = /3/2), which gives an 
estimate for the amount of spectral weight around the 
Fermi energy. This procedure was also used in Ref. 49. 
Additionally, the compressibility d (n) /dfi = Aq=o,tc;=o 
can be obtained directly from the susceptibility. Both 
quantities indicate that the normal phase is metallic. 

Appendix D: Approaching the transition 

In Fig. 7, we show the real-space correlation function 
as the ordered phases are approached. These figures are 
zoomed in, compared to the figures of the main text, only 
the sites close to the origin are shown. 
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